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We study the statistical properties of a simple genetic regulatory network that provides heterogeneity within 
a population of cells. This network consists of a binary genetic switch in which stochastic flipping between the 
two switch states is mediated by a "flipping" enzyme. Feedback between the switch state and the flipping rate is 
provided by a linear feedback mechanism: the flipping enzyme is only produced in the on switch state and the 
switching rate depends linearly on the copy number of the enzyme. This work generalises the model of [Phys. 
Rev. Lett. , 101, 1 1 8 104] to a broader class of linear feedback systems. We present a complete analytical solution 
for the steady-state statistics of the number of enzyme molecules in the on and off states, for the general case 
where the enzyme can mediate flipping in either direction. For this general case we also solve for the flip time 
distribution, making a connection to first passage and persistence problems in statistical physics. We show that 
the statistics of the model are non-Poissonian, leading to a peak in the flip time distribution. The occurrence of 
such a peak is analysed as a function of the parameter space. We present a new relation between the flip time 
distributions measured for two relevant choices of initial condition. We also introduce a new correlation measure 
to show that this model can exhibit long-lived temporal correlations, thus providing a primitive form of cellular 
memory. Motivated by DNA replication as well as by evolutionary mechanisms involving gene duplication, we 
study the case of two switches in the same cell. This results in correlations between the two switches; these can 
either positive or negative depending on the parameter regime. 

PACS numbers: 87.18.Cf,87.16.Yc,82.39.-k 



I. INTRODUCTION 

Populations of biological cells frequently show stochastic 
switching between alternative phenotypic states. This phe- 
nomenon is particularly well-studied in bacteria and bacterio- 
phages, where it is known as phase variation 1 1 1. Phase varia- 
tion often affects cell surface features, and its evolutionary ad- 
vantages are believed to involve evading attack from host de- 
fense systems (e.g. the immune system) and/or "bet-hedging" 
against sudden catastrophes which may wipe out a particu- 
lar phenotypic type. Switching between different phenotypic 
states is controlled by an underlying genetic regulatory net- 
work, which randomly flips between alternative patterns of 
gene expression. Several different types of genetic network 
are known to control phase variation — these include DNA 
inversion switches, DNA methylation switches and slipped 
strand mispairing mechanisms |[Tl|2][3l. 

In this paper, we study a simple model for a genetic net- 
work that allows switching between two alternative states of 
gene expression. Its key feature is that it includes a linear 
feedback mechanism between the switch state and the flipping 
rate. When the switch is active, an enzyme is produced and 
the rate of switching is linearly proportional to the copy num- 
ber of this enzyme. The statistical properties of this model 
are made non-trivial by this feedback, leading, among other 
things, to non-Poissonian behaviour that may be of advantage 
to cells in surviving in certain dynamical environments. Our 
model is very generic and does not aim to describe any spe- 
cific molecular mechanism in detail, but rather to determine 
in a general way the consequences of the linear feedback for 
the switching statistics. Motivated by the fact that cells of- 
ten contain multiple copies of a particular genetic regulatory 
element, due to DNA replication or DNA duplication events 



during evolution, we also consider the case of two identical 
switches in the same cell. We find that the two copies of the 
switch are coupled and may exhibit interesting and potentially 
important correlations or anti-correlations. Our model switch 
is fundamentally different from bistable gene networks that 
have been the subject of previous theoretical interest. In fact, 
as we shall show, our switch is not bistable but is intrinsically 
unstable in each of its two states. 

Before discussing our model in detail, we provide a 
brief overview of the basic biology of genetic networks and 
summarise some previously considered models for genetic 
switches. Genetic networks are interacting, many-component 
systems of genes, RNA and proteins, that control the func- 
tions of living cells. Genes are stretches of DNA (^1000 
base pairs long in bacteria), whose sequences encode partic- 
ular protein molecules. To produce a protein molecule, the 
enzyme complex RNA polymerase copies the gene sequence 
into a messenger RNA (mRNA) molecule. This is known as 
transcription. The mRNA is then translated (by a ribosome 
enzyme complex) into an amino acid chain which folds to 
form the functional protein molecule. The production of a 
specific set of proteins from their genes ultimately determines 
the phenotypic behaviour of the cell. Phenotypic behaviour 
can thus be controlled by turning genes on and off. Regula- 
tion of transcription (production of mRNA) is one important 
way of achieving this. Transcription is controlled by the bind- 
ing of proteins known as transcription factors to specific DNA 
sequences, known as operators, usually situated at the begin- 
ning of the gene sequence. These transcription factors may be 
activators (which enhance the transcription of the gene they 
regulate) or repressors (which repress transcription, often by 
preventing RNA polymerase binding). A given gene may en- 
code a transcription factor that regulates itself or other genes, 
leading to complex networks of transcriptional interactions 
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between genes. 

There has been much recent interest among both physical 
scientists and biologists in deconstructing complex genetic 
networks into modular units |4|, and in seeking to understand 
their statistical properties using theory and simulation HIS. 
Of particular interest is the fact that genetic networks are in- 
trinsically stochastic, due to the small numbers of molecules 
involved in gene expression |7, 8|. This can give rise to het- 
erogeneity in populations of genetically and environmentally 
identical cells f7|. For some genetic networks, this hetero- 
geneity is "all-or-nothing": the population splits into two dis- 
tinct sub-populations, with different states of gene expression. 
Such networks are known as bistable genetic switches: they 
have two possible long-time states, corresponding to alterna- 
tive phenotypic states. Well-known examples are the switch 
controlling the transition from the lysogenic to lytic states in 
bacteriophage A 19] (TO), and the lactose utilisation network of 
the bacterium Escherichia coli fllj. Several simple mecha- 
nisms for achieving bistability have been studied, including 
pairs of mutually repressing genes lfT2l[T3l . positive feedback 
loops (14] and mixed feedback loops (15). Such bistable ge- 
netic networks can allow long-lived and binary responses to 
short-lived signals — for example, when a cell is triggered by a 
transient signal to commit to a particular developmental path- 
way. 

Theoretical treatments of bistable genetic networks usually 
consider the dynamics of the copy number (or concentration) 
of the regulatory proteins involved. This affects the activa- 
tion state of the genes, which in turn influences the rate of 
protein production. The macroscopic rate equation approach 
lfT6l provides a deterministic (mean-field) description of the 
dynamics that ignores fluctuations in protein copy number or 
gene expression state. This approach, applied to a switch with 
two mutually repressing genes, has shown that co-operative 
binding of regulatory proteins is an important factor in gen- 
erating bistability [13|. Other studies have shown, however, 
that bistability can be achieved even when the deterministic 
equations have only one solution, due to stochasticity and 
fluctuations in protein numbers ifTTl [TSl . An alternative ap- 
proach is to study the dynamics of stochastic flipping between 
two stable states using stochastic simulations 1 19l l20ll2Tll . by 
numerically integrating the master equation |22|, or by path 
integral-type approaches |23|. This dynamical problem bears 
some resemblance to the Kramers problem of escape from a 
free energy minimum ll24l |25l . and one expects on general 
grounds that the typical time spent in one of the bistable states 
should be exponentially large in the typical number of pro- 
teins present in the state. This has been confirmed, at least 
for cooperative toggle switches formed of mutually repress- 
ing genes 1 19 20 1. From the perspective of statistical physics, 
interesting questions arise concerning the distribution of es- 
cape times and the connection to first passage properties of 
stochastic processes. 

In this paper, however, we are concerned with an intrinsi- 
cally different situation from these bistable genetic networks. 
The molecular mechanisms controlling microbial phase vari- 
ation typically involve a binary element that can be in either 
of two states. For example, this may be a short fragment of 



DNA that can be inserted into the chromosome in either of 
two orientations, a repeated DNA sequence that can be altered 
in its number of repeats, or a DNA sequence that can have two 
alternative patterns of methylation 1 1 1. The flipping of this el- 
ement between its two states is stochastic, with a flipping rate 
that is controlled by various regulatory proteins, the activity 
of which may be influenced by environmental factors. We 
shall consider the case where a feedback exists between the 
switch state and the flipping rate. This is particularly inter- 
esting from a statistical physics point of view because it leads 
to non-Poissonian switching behaviour, as we shall show. Our 
work has been motivated by several examples. The^m system 
in uropathogenic strains of the bacterium E. coli controls the 
production of Type 1 fimbriae (or pili), which are "hairs" on 
the surface of the bacterium. Individual cells switch stochas- 
tically between "on" and "off" states of fimbrial production 
lU |26l |22l EHl . The key feature of the fim switch is a short 
piece of DNA that can be inserted into the bacterial DNA in 
two possible orientations. Because this piece of DNA con- 
tains the operator sequence for the proteins that make up the 
fimbriae, in one orientation, the fimbrial genes are transcribed 
and fimbriae are produced (the "on" state) and in the other ori- 
entation, the fimbrial genes are not active and no fimbriae are 
produced (the "off" state). The inversion of this DNA element 
is mediated by recombinase enzymes. Feedback between the 
switch state and the switch flipping rate arises because the 
FimE recombinase (which flips the switch in the on to off 
direction), is produced more strongly in the on switch state 
than in the off state. This phenomenon is known as orienta- 
tional control |29, 30 31 1. The production of a second type of 
fimbriae in uropathogenic E. coli, Pap pili, also phase varies, 
and is controlled by a DNA methylation switch H] |2l [32l . 
Here, the operator region for the genes encoding the Pap pili 
can be in two states, in which the DNA is chemically modi- 
fied (methylated) at different sites, and different binding sites 
are occupied by the regulatory protein Lrp. Switching in this 
system is facilitated by the Papl protein, which helps Lrp to 
bind [33}. Feedback between the switch state and the flipping 
rate arises because the production of Papl itself is activated 
by the protein PapB, which is only produced in the "on" state 

A common feature of the above examples is the existence of 
a feedback mechanism: in the fim system this occurs through 
orientational control, and in the pap system, through activa- 
tion of the papl gene by PapB. In this paper, we aim to study 
the role of such feedback within a simple, generic model of 
a binary genetic switch. We shall assume that the feedback is 
linear, and we thus term our model a "linear feedback switch". 
In a recent publication |35|, we introduced a simple mathe- 
matical model of a DNA inversion genetic switch with orien- 
tational control, which was inspired by the fim system. Our 
model reduces to the dynamics of the number of molecules of 
a "flipping enzyme" R, which mediates switch flipping, along 
with a binary switch state. Enzyme R is produced only in the 
on switch state. As the copy number of R increases, the on 
to off flipping rate of the switch increases and this results in a 
non-Poissonian flipping process with a peak in the lifetime of 
the on state. The model is linear in the sense that the rate at 
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which the switch is turned off is a linear function of the num- 
ber of enzymes R which it produces, hi our previous work 
||35l . we imagined enzyme i? to be a DNA recombinase, and 
the two switch states to correspond to different DNA orien- 
tations, in analogy with the^m system. However, the same 
model could be used to describe a range of molecular mech- 
anisms for binary switch flipping with feedback between the 
switch state and flipping rate, and can thus be considered a 
generic model of a genetic switch with linear feedback. 

In our recent work 11351, we obtained exact analytical ex- 
pressions for the steady state enzyme copy number for our 
model switch with linear feedback, in the particular case 
where the flipping enzyme switches only in the on to off direc- 
tion (this being the relevant case for Jim). We also calculated 
the flip time distribution for this model analytically. Concep- 
tually, such a calculation is reminiscent of the study of persis- 
tence in statistical physics ll36l where, for example, one asks 
about the probability that a spin in an Ising system has not 
flipped up to some time [37 1. For the flip time distribution, 
we introduced different measurement ensembles according to 
whether one starts the time measurement from a flip event (the 
Switch Change Ensemble) or from a randomly selected time 
(the Steady State Ensemble). In the present paper, we extend 
this work to present the full solution of the general case of 
the model and extend our study of its persistence properties. 
The introduction of a rate for the enzyme mediated off to on 
flipping (kf^) has most significant effects on the flip time dis- 
tributions F{T), as illustrated in Figs. |6]and|7]where we show 
the parameter range over which a peak is found in F{T) for 
zero and non-zero kf^. We also prove an important relation 
between the two measurement ensembles defined in |35 | and 
use it to show that a peak in the flip time distribution only 
occurs in the Switch Change Ensemble and not in the Steady 
State Ensemble. We find that the non-Poissonian behaviour 
of this model switch leads to interesting two-time autocorre- 
lation functions. We also study the case where we have two 
copies of the switch in the same cell and find that these two 
copies may be correlated or anticorrelated, depending on the 
parameters of the model, with potentially interesting biologi- 
cal implications. 

The paper is structured as follows. In section 11 we de- 
fine the model, describe its phenomenology, and show that a 
"mean-field", deterministic version of the model has only one 
steady state solution. In section III we present the general so- 
lution for the steady state statistics and in section IV we study 
first passage-time properties of the switch; technical calcula- 
tions are left to the appendices. In section V we consider two 
coupled model switches and we present our conclusions in 
section VI. 



II. THE MODEL 

We consider a model system with a flipping enzyme R and 
a binary switch S, which can be either on or off (denoted re- 
spectively as 5*011 and S'otf)- Enzyme R is produced (at rate 
/c2) only when the switch is in the on state, and is degraded at 
a constant rate ki, regardless of the switch state. This repre- 
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FIG. 1 : (colour online) A schematic illustration of the model DNA 
inversion switch. 



sents protein removal from the cell by dilution on cell growth 
and division, as well as specific degradation pathways. Switch 
flipping is assumed to be a single step process, which can ei- 
ther be catalysed by enzyme R, with rate constants fc™ and 
fcg* and linear dependence on the number of molecules of R, 
or can happen "spontaneously", with rates /c™ and kf^. We 
imagine that the "spontaneous" switching process may in fact 
be catalysed by some other enzyme whose concentration re- 
mains constant and which is therefore not modelled explicitly 
here. Our model, which is shown schematically in Fig. [T] is 
defined by the following set of biochemical reactions: 
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A. Phenomenology 

We notice that there are two physically relevant and cou- 
pled timescales for our model switch: the timescale associ- 
ated with changes in the number of R molecules (dictated by 
the production and decay rates ki and k2), and that associated 
with the flipping of the switch (dictated by k^, k^ and the R 
concentration). 

We first consider the case where the timescale for R produc- 
tion/decay is much faster than the switch flipping timescale. 
The top left panel of Fig. |2] shows a typical dynamical tra- 
jectory for parameters in this regime. Here, we plot the num- 
ber n of R molecules, together with the switch state, against 
time. This result was obtained by stochastic simulation of re- 
action set ([T]i using the Gillespie algorithm I38. .39il . This al- 
gorithm generates a continuous time Markov process which 
is exactly described by the master equation ( 10 1. For a given 



switch state, the number n of molecules of R varies accord- 
ing to reactions (lai. When the switch is in the on state, n 
grows towards a plateau value, and when the switch is in the 
off state, n decreases exponentially towards n = 0. The time 
evolution of n can thus be seen as a sequence of relaxations 
towards two different asymptotic steady states, which depend 
on the switch position. To better understand this limiting case, 
we can make the assumption that the number of R molecules 
evolves deterministically for a given switch state. We can then 
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FIG. 2: (colour online) LEFT: Typical trajectories of the system when kf = kf^ — fca is increased (from top to bottom fcs = 0.0001, 0.01 
and 1). The other parameters are ki = 1, ^2 = 100 and fc™ = ^4'' = k4 — 0.1. Grey shading denotes periods in which the switch is in the 
on state, and the solid lines denote the number of enzyme molecules, plotted against time. In the bottom panel, the switch flips so fast that the 
grey shading is only shown in the inset where the trajectory from kit = [60, 61] is shown in detail. RIGHT: Probability distribution functions 
for the number nof R molecules, for parameter values corresponding to the trajectories shown in the left panels. The symbols are the result of 
numerical simulations (see text for details). The full curves plot the analytical results Eqs. \26\ and {36\ , which are in perfect agreement with 
the simulations. 



write down deterministic rate equations corresponding to the 
reaction scheme (jT]). These equations are first order differen- 
tial equations for p, the mean concentration of the enzyme. 
When the switch is on, the rate equation reads 



dp 



-kip + k2 



with solution 



p(0 = p(0)e-'=^* + g[l-e-'=^*] 



(2) 



(3) 



Thus the plateau density in the on state is given by the ratio 



(4) 



and the timescale for relaxation to this density is given by fci, 
the rate of degradation of Ri. When the switch is in the off 
state, the rate equation for p reads instead 



dp 
~dt 



-kip 



(5) 



and one simply has exponential decay to p = with de- 
cay time fci. In this parameter regime, switch flipping typ- 
ically happens when the number of molecules of R has al- 
ready reached the steady state (as in the top left panel of 
Fig. |2]l. Thus, the on to off switching timescale is given by 
l/(Pon^3" + ^4°); where pon is the plateau concentration of 
flipping enzyme when the switch is in the on state, given by 
Eq.(|4]). Since the corresponding plateau concentration in the 
off switch state is zero, the off to on switch flipping timescale 
is simply given by 1/ fc°*. 



We now consider the opposite scenario, in which switching 
occurs on a much shorter timescale than relaxation of the en- 
zyme copy number A typical trajectory for this case is shown 
in the bottom left panel of Fig. |2] Here, switching reactions 
dominate the dynamics of the model, and the dynamics of the 
enzyme copy number follows a standard birth-death process, 
with an effective birth rate given by the enzyme production 
rate in the on state multiplied by the fraction of time spent in 
the on state. A more quantitative account for these behaviours 



is provided later on, in III B 



For parameter values between these two extremes, where 
the timescales for switch flipping and enzyme number relax- 
ation are similar, it is more difficult to provide intuitive in- 
sights into the behaviour of the model. A typical trajectory 
for this case is given in the middle left panel of Fig. |2] Here, 
we have set the on to off and off to on switching rates to be 
identical: kf = kf^ and A;™ = /cf . We notice that typically, 
less time is spent in the on state than in the off state. As soon 
as the switch flips into the on state, the number of R molecules 
starts increasing and the on to off flip rate begins to increase. 
Consequently, the number of R molecules rarely reaches its 
plateau value before the switch flips back into the off state. 



To illustrate the effects of including the parameter kf^, 
we also show trajectories for different values of the ratio 
r — k'^^/kf in Fig. [s] for fixed fc™. For small r, the amount 
of enzyme decays to zero in the off state before the next off- 
to-on flipping event resulting in bursts of enzyme production. 
In contrast, when ?■ is 0(1), flipping is rapid in both directions 
so that p{n) is peaked at intermediate n. 
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FIG. 3: (colour online) LEFT: Typical trajectories of the system when r = kf/kf^ = is increased (from top to bottom r = 0, 0.5 and 1). The 
other parameters are ki — 1, k2 = 100, fc™ = 1 and fc™ = k^^ = ^4 = 0.1. Grey shading denotes periods in which the switch is in the on 
state, and the solid lines denote the number of enzyme molecules, plotted against time. In the bottom panel, the switch flips so fast that the 
grey shading is only shown in the inset where the trajectory from kit = [60, 61] is shown in detail. RIGHT: Probability distribution functions 
for the number nof R molecules, for parameter values corresponding to the trajectories shown in the left panels. The symbols are the result of 
numerical simulations (see text for details). The full curves plot the analytical results Eqs. \26\ and {36\ , which are in perfect agreement with 
the simulations. 



B. Mean-field equations 

To explore how the switching behaviour of our model 
arises, we can write down mean-field, deterministic rate equa- 
tions corresponding to the full reaction scheme Q. These 
equations describe the time evolution of the mean concen- 
tration p{t) of R molecules and the probabilities Qon{t) and 
Qofi{t) of the switch being in the on and off states. These 
equations implicitly assume that the mean enzyme concentra- 
tion p is completely decoupled from the state of the switch. 
Thus correlations between the concentration p and the switch 
state are ignored and the equations furnish a mean-field ap- 
proximation for the switch. As we now show, this crude 
type of mean-field description is insufficient to describe the 
stochastic dynamics of the switch, except in the limit of high 
flipping rate. Noting that (Jon (i) + (? off (^) — 1, the mean-field 
equations read: 



dpjt) 
dt 



k2Qon{t) - kip(t) 



(6a) 



dQonjt) 

dt 



= {kf + p{t)kf){l - Q,„it)) 

- {kT + pit)k"^nQonit) 



(6b) 



The above equations have two sets of possible solutions for the 
steady-state values of p and Qon, but only one has a positive 
value of p, and is therefore physically meaningful. The result 
is: 



n k' 



.off 



(fcf + fcf ) + VA 



2{kf + /C3™) 



(7) 



where 

A = (ponfcf - (fcf + fc^))" + 4po„fcf (fcf + kT) 



and 



Qon = pI Po 



(8) 



(9) 



The most interesting conclusion to be drawn from this mean- 
field analysis is that there is only one physically meaningful 
solution. In this solution, the enzyme concentration p is less 
than the plateau value in the on state [pon of Eq.(|4]i]. Thus 
reaction scheme ([TJ does not have an underlying bistability. 
The two states of our stochastic switch evident in Figures [2] 
and|4]for low values of k^ and fc4 are not bistable states but are 
rather intrinsically unstable and transient states, each of which 
will inevitably give rise to the other after a certain (stochasti- 
cally determined) period of time. In this sense, our model is 
fundamentally different from the bistable reaction networks 
which have previously been discussed lT3l [191 1401 . On the 
other hand, in the limit of rapid switch flipping, where fca or 
fc4 is large, the mean-field description holds and the protein 
number distribution does show a single peak whose position 
is well approximated by Eq. (j?]), as shown in Figures |2] and |4] 
for the case k^ = 1. 



III. STEADY-STATE STATISTICS 

A. Analytical solution 

Returning to the fully stochastic version of the reaction 
scheme ([T]i, we now present an exact solution for the steady- 
state statistics of this model. A solution for the case where 
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FIG. 4: (colour online) LEFT: Typical trajectories of the system when fc™ — kf' — is increased (from top to bottom ki — 0.1, 1 and 
100. Other parameters are fci = 1, k2 = 100 and fc™ — k°f ~ ks — 0.001. In each panel the grey shading denotes that the switch is on 
and the line plots the number of enzymes against time. In the third panel the grey shading is only shown in the inset where the trajectory 
from kit — [60, 61] is detailed. RIGHT: Probability distribution functions of the number of Ri molecules in the cell for parameter values 
corresponding to the trajectories shown in the left panels. The symbols are the result of numerical simulations (see text for details). The full 
curves plot the analytical results Eqs. \26\ and \36\ and pass perfectly through the simulation points. 



fcg* = was sketched in Ref. f35\. Here we present a com- 
plete solution for the general case where kf^ ^ 0, and we 
discuss the properties of the steady-state as a function of all 
the parameters of the system. 

We first define the probability Ps(n, t) that the system has 
exactly n enzyme molecules at time t and the switch is in the 
s state (where s = {on, off}). The time evolution of ps is 
described by the following master equation: 



dpsin) 



{n+l)kiPs{n+l)+k^Ps{n-l)+nkl ''pi_s(n) 



dt 

+ kl^'pi^sin) - {nki + k^ + nk^ + kl)ps{n) , (10) 

where we use the shorthand notations {off, on} = {0,1}, 
/c2* = and /c™ = k2- In the steady state, the time deriva- 



tive in Eq.( 10 1 vanishes, and the problem reduces to a pair of 
coupled equations for pon and poff. 

{n+l)kiPo»{n+l)+k2Pon{n-l)+nk°^Pos{n)+kf^Pofi{n) 
= (nfci + fc2 + nkf + k°^)pM , (11a) 

(n + l)/ciPoff(" + 1) + nk°^p„^{n, t) + fc°>on(f^, t) 

^ {nki + nkf + fcf )poff(n, t) . (lib) 

To solve the above equations we introduce the generating 
functions 



Gs{.z) 



oo 



n=0 



Ps{n)z''' 



(12) 



The steady-state equations (111 can be now written as a set of 



linear coupled differential equations for Gg'- 

CiGon{z) = C2Goff{z) , 
Cs,Gofs{z) = CiGoniz) , 



(13a) 
(13b) 



where Ci are linear differential operators: 

Li{z) =ki{z - l)d, ~ k2{z - 1) + A:™z9, + kf , (14a) 

C2{z) ^kfzd, + kf , (14b) 

C:,{z) =ki{z - 1)9, + fcf za, + kf , (14c) 

Ci{z) ^k°^zd, + k°l . (14d) 



In order to solve the two coupled Eqs. ( [T3] l it is first useful 
to take their difference. After simplification this yields the 
relation: 



dzGofi{z) = -dzGon{z) + -j-Gon{z) 

fcl 



(15) 



Next, we take the first derivative of (TTSbji and then replace the 
derivatives of Goff with the relation ( 15 1. After some algebra. 



one finds that Gon verifies the following second order differ- 
ential equation: 

ki{az~ki)G':„{z) + {kiP—fz)G'^^{z)-5G,,{z) = , 

(16) 

where the Greek letters are combinations of the parameters of 
the model: 



a = ki + fc™ 

P=ki+k2 

7 = fc2(fci + k 



hfcf 
kf, 

offN 

3 ) ; 
I off 



roff 



S = fc2(fci + kf + kf) . 
We now introduce the new variable 



u{z) 



1 

fcl a 



1_ 
a2 



UO + z{ui - Mo) 



(17a) 
(17b) 
(17c) 
(17d) 

(18) 
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and the new parameter combinations: 

C = ^^o+-, V=- ■ (19) 
a 7 

We can now write Gon(z) (and Gos{z)) in terms of the vari- 



(20) 



able w (1 8 1 by defining the functions 

Js{u) = G,(z) . 
The differential equation ^TE[ then reads: 

+ (C - u)JM - vJM = . (21) 

Looking for a regular power series solution of the form 



Jon{u) = ^ a„u" 



n=0 



one obtains the following solution: 

Jon(u) = flo iFi (r/, C,u) , 



(22) 



(23) 



where li^i denotes the confluent hypergeometric function of 
the first kind, 



(On n\ 



(24) 



and {a)n — a{a+\) . . . (a + n— 1) denotes the Pochhammer 
symbol. 

The constant ao will be determined using the boundary con- 
ditions, which we discuss later We first note that the above 
result for Jon{u) can be translated into Gon{z) by replacing u 
with the expression of u{z) in (22 1 and expanding in powers 
of z: 



Gon{z) = ^ a„(wo + z{ui - uo))" 

n=0 
oo n 

n—0 7n—0 
oo oo 



n—0 m—n 



(25) 



where we have relabelled the indices n — m ^ n and n 



m in the last line. We can identify Ponin) from (12 1 as the 
coefficient of z" in the above expression: 



nj = 2^ OmUg (ui - uo) 

m—n 

From ( |22] i and ( |23] l we read off 

_ Oo (??)„ 



(26) 



(27) 



Substituting ( 27 1 in ( 26 1 we deduce, using the definition of the 

hypergeometric function (24 1 and noting (Q;)„+m ~ 
n)m, that 

= ao j ——iFi{r] + n,C + n,uo). (28) 

J^! (On 



In deriving this expression we have, in fact, established the 
following identity which will prove useful again later: 



n=0 



- Up)'" {ri)n 
n\ (On 



(29) 

To compute Goff(2:), we integrate Eq.( 15 1, which yields, using 
the form of Jon(u) (23 1: 



Goii{z) + Gon(2) 
fc2(C-l) 



ao 



fci(77 - l)(ui - Uo) 



li^i (77 -1,C-1, «.) = «:, (30) 



where k is our second integration constant. We then have 
two constants, ao and k, which still need to be determined. 
The constant k can be found using the normalisation con- 
dition 'YlirXVoxiip) + Poff("-)) = 1. which is equivalent to 
Gon(l) + Goff(l) = 1. Using this condition, we obtain 

^ = l-ao , , -ii^i(r;-l,C-l,iii) . (31) 

fci(77 - l)(wi - Uo) 

In order to compute the remaining constant ao, we consider 
the boundary condition at z = 0. From the definition ( [T2| of 
the generating function we see that Gs{z = 0) = -psin = 0). 
Our boundary condition thus reads: 



Jon(uo) + Joff(wo) = Pon(O) + Poff (0) 



(32) 



Setting n = in the master equation Eq.( |l la| l [noting that the 
term in Pon{n — 1) vanishes] gives Poff(O) in terms of Pon(O) 
and Pon(l): 



Poff(O) = 



-Pon(O) - 7TfPon(l) . 



(33) 



Combining Eqs.(30l [with z — 0] and ( [3T] l, substituting in 
Eq.(32 1, using Eq.(33 i to eliminate Poff(O), and finally substi- 
tuting in expressions for Pon(O) and Pon(l) from Eq.(26l, we 
determine ao: 



On ^ = f 1 + ■ 1 li^l(?7, (,""0) 



,off 



kiVjui - Up) 

kfc 



Fi(?7+l,C + l,uo) 



[iFi(7/-l,C-l,uo) 



ki{ri - 1)(mi - Uo) 

- iFi(77-l,C-l,?/i)] . (34) 

The final step in obtaining our exact solution is to provide an 
explicit expression forpoff(n). From (30i we have 



Goff(z) = K - Gon{z) 

k2{C-l) 



ap 



fci(?7 - l)(ui - Uo) 



iF,{7^-lX-l,u,) , (35) 
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and using the identity (j29]l we obtain: 

Poff{n) = KSnfi + 



{ui - uoY 



(Or 



-n-lX 



n- l,Uo) 
(36) 



where di , is the Kronecker delta. 



Our exact analytical solution ( 26 1, ( 34 1 and ( 36 1 is verified 
by comparison to computer simulation results in the right pan- 
els of Figs. |2]and|4] Here, we plot the probability distribution 
function for the total number of enzyme molecules; 



p{n) = Ponin) +Pofi{n) . 



(37) 



Computer simulations of the reaction set ^ were carried 
out using Gillespie's stochastic simulation algorithm |38, 39|. 
Perfect agreement is obtained between the numerical and an- 
alytical solutions, as shown in Figs. |2]and|4] 



B. Properties of the steady-state 

Having derived the steady-state solution forp(n), we now 
analyse its properties as a function of the parameters of the 
model. We choose to fix our units of time by setting ki, the 
decay rate of enzyme R, to be equal to unity (so our time 
units are fcf ^). With these units, the plateau value for the 
number of enzyme molecules in the on switch state is given 
by Pon — In this section, we will only analyse the case 
where pon — 100. To further simplify our analysis, we set 



Loff 

uoff 



/c.s and feS 



Loff 



where kf = and fc°' 



ki (a discussion of the case 
7^ is provided in Ref. 1 35 1). We then 
analyse the probability distribution p{n) as a function of the 
i?-dependent switching rate k^ and the i?-independent switch- 
ing rate fc4. The results are shown in the right-hand panels of 
Fig. |2]and Fig. |4] We consider the three regimes discussed 
in section II A| that in which enzyme number fluctuations are 



much faster than switch flipping, that where the opposite is 
true, and finally the regime where the two timescales are sim- 
ilar. 

In the regime where switch flipping is much slower than en- 
zyme production/decay [fci ^ (fc™ + k2k™ /ki)}, the proba- 
bility distribution p{n) is bimodal. This is easily understand- 
able in the context of the typical trajectories shown in the left 
top panels in Figs. |2] and |4] in this regime, the number of 
molecules of R always reaches its steady-state value before 
the next switch flip occurs. It follows then that Pon{n) is a 
bell-shaped distribution peaked around k2/ki, while Posin) 
is highly peaked around zero, so that the total distribution 
p{n) ^ Pon{n) + Pos{n) is bimodal. 

In contrast, when switching occurs much faster than en- 
zyme number fluctuations the probability distribution p{n) 
is unimodal and bell shaped, as might be expected from the 
trajectories in the bottom left panels of Figs. |2]and|4] As 
discussed in section II A in this regime the number of R 



molecules behaves as a standard birth-death process with ef- 
fective birth rate given by k2 multiplied by the average time 
the switch spends in the on state, and death rate ki . For such a 
birth-death process the steady state probability p{n) is a Pois- 
son distribution with mean given by the ratio of the birth rate 
to the death rate. To show that our analytical result reduces to 
this Poisson distribution, we consider the case where enzyme- 
mediated switching dominates (as in Fig. |2ji, so that both kf^ 
and fcg" are much greater than fci. The fraction of time spent 
in the on state is k^^/ (/c™ + kf^) , thus the effective birth rate 
is k2kf/ (kf' + kf). In the Umit A:§" -> oo and kf -> oo 
with r = kf/k'^" constant, one finds that rj ^ 1, ( 1, and 
Uz k2rz/[ki{l+r)]. Using the fact that ii^i(l,l,a;) — e^, 
Eq.([23]l gives, in this limit. 



Gon(z) = aoexp 



k2rz 
fci(l + r) 



(38) 



which is the generating function of a Poisson distribution 
with mean k2kf/[ki{k™ + kf)]. Plugging this result into 
Eq.((30| and taking again the limit k^ ^ oo [and using 
that iFi{0,0,x) = 1] finally yields the result that p{n) = 
Pon{n) + Posin) is indeed a Poisson distribution. The same 
approach can be taken for the case of Fig. |4] where k^ is 
constant, and fc™ and kf become very large. The probabil- 
ity distribution p{n) then becomes a Poisson distribution with 
mean k2kf /[ki{k1" + kf)]. The above result is only valid 
when r 7^ 0. In fact, as shown in Fig. [3] when r = the dis- 
tribution of R is peaked at and does not have a Poisson-like 
shape. 

Finally, when there is no clear separation of timescales be- 
tween enzyme number fluctuations and switch flipping, the 
distribution function for the number of enzyme molecules has 
a highly non-trivial shape, as shown in the middle panels of 
Figs.|2]and|4] 



IV. FIRST PASSAGE TIME DISTRIBUTION 

We now calculate the first passage time distribution for our 
model switch. We define this to be the distribution function 
for the amount of time that the switch spends in the on or off 
states before switching. This distribution is biologically rele- 
vant, since it may be advantageous for a cell to spend enough 
time in the on state to synthesise and assemble the components 
of the "on" phenotype (for example, fimbriae), but not long 
enough to activate the host immune system, which recognises 
these components. The calculation for the case kf = was 
sketched in |351. Here we provide a detailed calculation of the 
flip time distribution in the more general case ^53'^ 7^ 0. We 
find that this dramatically reduces the parameter range over 
which the flip time distribution has a peak. We demonstrate 
an important relation between the flip time distributions for 
the two relevant choices of initial conditions (Switch Change 
Ensemble and Steady State Ensemble). The first passage 
time distribution is important and interesting from a statistical 
physics point of view as it is related to "persistence". Gener- 
ally, persistence is expressed as the probability that the local 
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value of a fluctuating field does not change sign up to time 
t ll36l . For the particular case of an Ising model, persistence 
is the probability that a given spin does not flip up to time t. 
In our model, the switch state S plays the role of the Ising 
spin. For other problems, there has been much interest in the 
long-time behaviour of the persistence probability, which can 
often exhibit a power-law tail. In our case, however, we ex- 
pect an exponential tail for the distribution of time spent in 
the on state, because linear feedback will cause the switch to 
flip back to the off state after some characteristic time. We 
are therefore interested not only in the tail of the first passage 
time distribution, but in its shape over the whole time range. 

A. Analytical results 



To compute the distribution F{T), we first consider the sur- 
vival probability (n, t), that, given that at time t = (cho- 
sen according to ensemble W), the switch was in state s, at 
time t it is still in state s and has n molecules of enzyme R. 
As the ensemble W only enters through the initial condition, 
we may drop the superscript W in what follows. The evolu- 
tion equation for hg is the same as for ps{n,t), but without 
the terms denoting switch flipping into the s state. This re- 
moves the coupling between po„ and poff that was present in 
the evolution equations ([TT|)): 



Tr:Kn{n,t) = (n + l)kihon{n + l,t) + k2Kn{n - l,t) 
at 

- {nki + k2 + nkf + k°^)K,,{n,t) , (42a) 



We consider the probability Fs{T\nQ)dT that if we begin 
monitoring the switch at time to when there are molecules 
of the flipping enzyme R, it remains from time ^ ^ + T 
in state s, and subsequently flips in the time interval t^ + T ^ 
^0 + T + dT. This probability is averaged over a given en- 
semble of initial conditions, determined by the experimental 
protocol for monitoring the switch. Mathematically, the initial 
condition rtg for switch state s is selected according to some 
probability VFs(no) and we define 



F,(T) =^F,(T|no)W^.(no) 



as the flip time distribution for the ensemble of initial condi- 
tions given by ^^^ (no). 

The most obvious protocol would be to measure the interval 
T from the moment of switch flipping, so that the times to 
correspond to switch flips and the T are the durations of the on 
or off switch states. We call this the Switch Change Ensemble 
(SCE). In this ensemble, the probability 14^^*^^ of having n 
molecules of R at the time <o when the switch flips into the s 
state is: 



pi-s{n){nkl ' + kl 



l-s 
3 



K-1 



(40) 



where for notational simplicity, s = {I7O} represents 
{on, off}. The numerator of the r.h.s of Eq.(40l gives the 



steady state probability that there are n molecules present in 
state 1 — s, multiplied by the flip rate into state s. The denom- 
inator normalises Wg'^^{n). 

We also consider a second choice of initial condition, which 
we denote the Steady State Ensemble (SSE). Here, the initial 
time io is chosen at random for a cell that is in the s state. 
This choice is motivated by practical considerations: experi- 
mentally, it is much easier to pick a cell which is in the s state 
and to measure the time until it flips out of the s state, than 
to measure the entire length of time a single cell spends in the 
s state. The probability W^^^ of having n molecules of R at 
time tQ is then the (normalised) steady-state distribution for 
the s state: 



W. 



SSE 



Ps{n) 



(41) 



d_ 
di 



hofi{n, t) = {n + l)fci/ioff(n + 1, 



{nki + nkf + kf)Kg{n,t) . (42b) 
Introducing the generating function 



hs{z, = ^ z''hs{n, t) 



(43) 



n=0 



Qc^-^ the above equations reduce to: 



d ~ 

^^K,{z,t) ^ {ki - {ki + kf)z)d,K,{z,t) 

+ {k2Z-{k2 + k°^'))h,,{z,t) , (44a) 



^h,s{z.t) = (/ci - (fci + kf)z)dX«{z,t) 

-kfhos{z,t) . (44b) 

We can relate to F by noting that J2n hs{n, t) = hs{l,t) 
is the total probability that the switch has not flipped up to 
time t. Hence, 



(45) 



Equations (44 1 can be solved using the method of characteris- 
tics [,41 j . The result, detailed in Appendix [A] is: 

X W?(A:iTon + e-*/^"(z-fciro„)) , (46) 



and Uon = A:™ + ^2(1 - kiTon)- 



where Ton = (fci + fc™ 
The function W is the generating function for the distribution 
of enzyme numbers W{n) at the starting time for the mea- 
surement: 



(47) 



where W refers to W^'"^ or W^^^. The function hof{{z,t) 
can be obtained in an analogous way: this produces the same 
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expression as for hon, but with k2 set to zero and with all "on" 
superscripts replaced by "off: 

Kff{z, t) = e-<*VK(fciroff + e-'/^-^Cz - hr^s)) , (48) 

so that Toff = (fci + kf^)^^. We can then obtain the distribu- 
tions Fo„{T) and Foff(r) by differentiating the above expres- 
sions, according to Eq.([45]i: 



ki]W' kiT, 



-T/r„, 



(1 - /ciTo, 




(49) 



Foff(T) = exp (-(fcf + l/roff)r) 



(50) 



In the above expressions, the function Ws is given for the 
steady state ensemble (SSE) by 



Wf^ = G,(z)/G,(l) 
and for the switch change ensemble (SCE) by 



(51) 



wSCE/ ^ _ ^3 ''zG[_^{z) + kl ''G'i-s(z) 



B. Relation between SSE and SCE 

We now show that a useful and simple relation can be de- 
rived between Fsse{T) and i^scE(2^)- Let us imagine that we 
pick a random time t, chosen uniformly from the total time 
that the system spends in state s. The time t will fall into an 
interval of duration T, as illustrated in Fig. |5] We can then 
split the interval T into the time Ti before t and the time T2 
after t, such that Ti + Ts = T. 

We first note that the probability that our randomly chosen 
time t falls into an interval of length T is: 



Prob(r) dT = 



TF^'-^{T) dT 



(53) 



Eq.(53 1 expresses the fact that the probability distribution for 
a randomly chosen flip time T is F^'^^{T) dT, but the prob- 
ability that our random time t falls into a given segment is 
proportional to the length of that segment. Since the time T 



off 



switch position 



Ti ' T2 
< ^ > 



t 



time 



Fo„(T) = exp (-(Won + 1/to„)T + fc2To„(l - e-^/"™)) 



FIG. 5: Schematic illustration of a possible time trajectory for the 
switch; t is a random time falling in an interval of total length T and 
splitting it into two other intervals denoted Ti and T2, as discussed 
in Section llVBl 



is chosen uniformly, the probability distribution for T2, for a 
given T, wiU also be uniform (but must be less than T): 



Prob(r2|r) dT - 



e(T - T2) 

T 



dT 



(54) 



One can now obtain F^ from Prob(r2|T) by integrating 



Eq.(54i over all possible values of T, weighted by the rela- 
tion (53 I. This leads to the following relation between F^ce 
andF^: 



J^^ Ff^^jT') dT' 



(55) 



Taking the derivative with respect to T2 this can be recast as 



dT 



(T) 



(56) 



SCE 



where (T)gcE is simply the mean duration of a period in the 
on state. We have verified numerically that the expressions 
(49)) and (|50l) for Ff^^{T) and Ff^^{T) derived above do in- 

This relation can also be under- 



deed obey the relation ( [56] l 
stood in terms of backward evolution equations as we discuss 
in Appendix |B] 



C. Presence of a peak in F{T) 

We now focus on the shape of the flip time distribution 
F{T), in particular, whether it has a peak. A peak in ^^^^^(r) 
could be biologically advantageous for two complementary 
reasons. Firstly, after the switch enters the on state there may 
be some start-up period before the phenotypic characteristics 
of the on state are established, so it would be wasteful for 
flipping to occur before the on state of the switch has become 
effective. Secondly, the on state of the switch may elicit a neg- 
ative environmental response, such as activation of the host 
immune system, so that it might be advantageous to avoid 
spending too long a time in the on state. For example, in the 
case of the^m switch, a certain amount of time and energy is 
required to synthesise fimbiiae, and this effort will be wasted 
if the switch flips back into the off state before fimbrial syn- 
thesis is complete. On the other hand, too large a population 
of fimbriated cells would trigger an immune response from the 
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host, therefore the length of time each cell is in the fimbriated 
state needs to be tightly controlled. We note that for bistable 
genetic switches and many other rare event processes, wait- 
ing time distributions are exponential (on a suitably coarse- 
grained timescale). This arises from the fact that the alter- 
native stable states are time invariant in such systems. The 
presence of a peak in F^^^{T) for our model switch would 
indicate fundamentally different behaviour, which occurs be- 
cause the two switch states in our model are time-dependent. 

The presence of a peak in the distribution F{T) requires 
the slope of F{T) at the origin to be positive. Applying this 
condition to the function Fo„ ( 49 1 we get: 



0.5 



- ik°''fW"{l) > . (57) 



Eq.(47i allows us to expressing the derivatives of W{1) as 



functions of the moments of n, so that we finally get our con- 
dition as a relation between the mean and the variance of the 
initial ensemble: 



where 
of Eq. 



(fc4 



A:™(fci + 2kl 



'Won 



- (fc™)' > ' (58) 

denotes an average taken using the weight Won 



(40i or (41 



Analogous conditions can be found for 
a peak in the off to on waiting time distribution. The mo- 
ments involved in the above inequality can be computed using 
the exact results of the previous section. The l.h.s. of (58 i 
can then be computed numerically for different values of the 
parameters, to determine whether or not a peak is present in 

Fin 

For the SSE, there is never a peak in the flip time distribu- 



tion. This follows directly from the relation ( 56 1 between the 
SSE and SCE, which shows that the slope of Ff^{T) at the 
origin is always negative: 



dT 



T=0 



f!^ 



< 



Thus a peak in the waiting time distribution cannot occur 
when the initial condition is sampled in the steady state en- 
semble. 

For the SCE, we tested inequality ( [58] l numerically and 
found that a peak in the distribution F{T) is possible for the 
time spent in the on state (F^jf^), but not for the off to on 
waiting time distribution (F^^^). This is as expected and can 
be explained by noting that to produce a peak in F^'^^{T), 
the flipping rate must increase with time in state s. In the on 
state the flipping rate typically does increase with time as the 
enzyme R is produced, while in the off state the flipping rate 
decreases in time as R decays. 

We now discuss the general conditions for the o ccurrence 
of a peak in F^^'^^ 



lllB 



that in the 



We first recall from section ] 
regime where the copy number of the enzyme R relaxes much 
faster than the switch flips [ki ^ fc™ + k2k™/ ki], the plateau 
level of R is reached rapidly after entering the on state, so 




FIG. 6: Occurrence of a peak in the waiting time distribution sam- 
pled in the Switch Change Ensemble. The shaded area delimits the 
region where there is a peak (here the parameters are: ki = 1, 
k2 = 10 and kf = fc^" = ks and kf = kf = k^). The dashed 
line dehmits the same region for k2 = 100. The insets show an in- 
stance of the distribution both in the SCE (solid red line) and in the 
SSE (blue dashed line): (a) There is a peak (fe = 10, fea — 0.1, 
k4 — 0.1); (b) On the transition line, where the slope at the origin 
vanishes (k2 = 10,A:3 = 0.15, ^4 = 0.209384...); (c) There is no 
peak (k2 = 10, fca = 0.2, k4 = 0.35). 



that the flipping rate out of the on state is essentially constant. 
This leads to effectively exponentially distributed flip times 
from the on state, so that no peak is expected. In the opposite 
regime, where switch flipping is much faster than R number 
relaxation [k^ 3> 0], we again expect Poissonian statistics and 
therefore exponentially distributed flip times. Thus it will be 
in the intermediate range of k^ that a peak in the flip time dis- 
tribution may occur The exact condition for this (fSSj) is not 



(59) particularly transparent as the dependence on the parameters 



is implicit in the values of the (n)^- and (?^^)^ . In particu- 
lar, the effects of the parameters k^ and fc2 are coupled, since 
the effective i?-mediated switching rate depends on the copy 
number of R. However we can make a broadbrush descrip- 
tion of what is required. First the switch should enter the on 
state with typical values of n <C Pon so that there is an initial 
rise in the value of n and therefore the flipping rate. Second, 
we expect that the flipping should be predominantly effected 
by the enzyme R rather than spontaneously flipping i.e. k^ 
should govern the flipping rather than k4. 

Fig. shows the region in the ki,-ki plane where F^^^ 
has a peak, for the case where fc™' = A;§" = /ca and fc™ = 
A;4. These results are obtained numerically, using the 



uoff 



inequality (|58|l. The distribution F^^^ is peaked for parameter 



values inside the shaded region. The insets show examples of 
the distributions F^^^{T) and F^^^{T) for vai'ious parameter 
values. At the boundary in parameter space between peaked 
and monotonic distributions (solid line in Fig. |6]l, ^^^'^^(r) 
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FIG. 7: Same plot as Fig. [6] but for ^3'^ — 0. The shaded area 
delimits the values of ^4 and fc™ (with ^2 = 10) for which there is 
a peak in the flip time distribution. The dashed line is the separation 
line for k2 = 100. The examples in the insets have as parameters 
k2 = 10 and: (a) k"^" = 15, ^4 = 0.15; (b) kf = 50, ^4 = 
0.162383...; (c) kf = 80, fc4 = 0.4. 



has zero gradient at T = (inset (b)). The dashed line in Fig. 
|6|l shows the position of the boundary for a larger value of 
the enzyme production rate k2- As ^2 increases, the range of 
values of for which there is a peak decreases. Increasing ^2 
increases the number of enzyme present, which will increase 
both the off to on and on to off switching frequency, since here 
^3" ^ ^3'^' ^ ^3- Thus it appears that approximately the same 
qualitative behaviour can be obtained for smaller values of 
when k2 is increased. 

In our previous paper f35|, we analysed the case where 
kf^ = 0: i.e. the flipping enzyme R switches only in the 
on to off direction. This case applies to the^w system. Fig. 
[Tjshows the analogous plot, as a function of fc™ and fc4, when 
^3*' = 0. The region of parameter space where a peak occurs 
in F^^^^{T) is much wider than for nonzero kf^. In this case 
an increase of k2 produces a larger range of parameter values 
A:™ for which there is a peak (dotted line in Fig. [Tjl. Here, the 
off to on switching process is i?-independent, and is mediated 
by /c4 only (since kf^ = 0). The typical initial amount of R 
present on entering the on state is thus not much affected by 
/c2, although the plateau level of R increases with k2. There- 
fore, as fc2 increases, the enzyme copy number in the on state 
becomes more time-dependent, increasing the likelihood of 
finding a peak. 

The comparison between Figs. |6]and|7]suggests that the rel- 
ative magnitudes of the i?-mediated switching rates in the on 
to off and off to on directions, fc™ and k^^, play a major role in 
determining the parameter range over which F^^^ is peaked. 
This observation is confirmed in Fig. |8] where the boundary 
between peaked and unpeaked distributions is plotted in the 
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FIG. 8: Diagram showing the occurrence of a peak when the ratio 
r — kf'/ fc™ is varied. Here ki = 1 and k2 = 10. The inset shows a 
zoom of the plot in the vicinity of fc™ = 0. 



kf'-k4 plane for various ratios r = kf^/k'^\ The larger the 
ratio r, the smaller the region in parameter space where there 
is a peak. An intuitive explanation for this might be that as r 
increases, the the typical initial number of R molecules in the 
on state increases, so that less time is needed for the R level to 
reach a steady state, resulting in a weaker time-dependence of 
the on to off flipping rate and less likelihood of a peak occur- 
ring in F{T). If the presence of a peak in F^^^ is indeed an 
important requirement for such a switch in a biological con- 
text, then we would expect that a low value of kf^, as is in fact 
observed for the^m system, would be advantageous. 



V. CORRELATIONS 

A peaked distribution of waiting times is by no means the 
only potentially useful characteristic of this type of switch. In 
this section, we investigate two other types of behaviour that 
may have important biological consequences: correlations be- 
tween successive flips of a single switch, and correlated flips 
of multiple switches in the same cell. We analyse these novel 
phenomena using numerical methods. We introduce a new 
correlation measure which enables us to quantify the extent 
of the correlation as a function of the parameter space. Our 
main findings are that a single switch shows time correlations 
which appear to decay exponentially, and that two switches in 
the same cell can show correlated or anti-correlated flipping 
behaviour depending on the values of k^^ and fc™. 
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A. Correlated flips for a single switch 

Biological cells often experience sequences of environmen- 
tal changes: for example, as a bacterium passes through the 
human digestive system it will experience a series of changes 
in acidity and temperature. It is easy to imagine that evolution 
might select for gene regulatory networks with the potential 
to "remember" sequences of events. The simple model switch 
presented here can perform this task, in a very primitive way, 
because it produces correlated sequences of switch flips: the 
amount of R enzyme present at the start of a particular period 
in state s depends on the recent history of the system. In con- 
trast, for bistable gene regulatory networks, or other bistable 
systems, successive flipping events are uncorrected, as long 
as the system has enough time to relax to its steady state be- 
tween flips. 

In our recent work ll35l . we demonstrated that successive 
switch flips can be correlated for our model switch, and that 
this correlation depends on the parameter k'^^: correlation in- 
creases as kf^ increases. Here, we extend our study and in- 
troduce a new measure of these correlations: the two time 
probability p{s, t; s', t') that the switch is in position s at time 
t and in position s' at time t'. In the steady state the two-time 
probability depends only on the time difference t = t — t'. In 
order to compare different simulations results, we define the 
auto-correlation function: 



C(r) 



Pon— on 



1, 



Poff 



(60) 



where Pon-oa{r) = p{on,t;on,t + t), _Poff-off(T") = 
p(off, t; off, t + t), and pon (Poff) is the probability of being in 
the on (off) state. The correlation function ( [60) 1 takes values 
between —1 and 1, in such a way that it is positive for positive 
correlations, negative for negative correlations and vanishes if 
the system is uncorrected. This function allows us to under- 
stand whether, given that the switch is in a given position s at 
time t, it will be in the same state s at a later time t + t. 
Fig. l9] shows simulation results for different values of 
= fcg = fca and = ^4*^ = fc4. As expected, the cor- 
relation function vanishes in the limit of large t, meaning that 
in this limit there are no correlations. Furthermore, we can 
see that the strength of the correlations decreases when either 
fcs or fc4 are increased. This is consistent with the previous 
remark that in the limit of large switching rate (i.e. either k^ 
or fc4) the distribution of enzyme numbers tends to a Poisson 
distribution. It is thus not surprising that in this same limit the 
correlations vanish. In the insets of Fig. |9]we plot the same 
correlation function on a semi-logarithmic scale. The data for 
the highest values of k^ or fc4 (the dotted green curves) is not 
shown since the decrease is too sharp, and does not allow for a 
clear interpretation. For the smallest values of ^3 and k4 (blue 
curves), the decay seems to be exponential. However, for in- 
termediate values of ^3 or k^ (dashed red curves) the evidence 
for an exponential decay is less clear and the issue deserves a 
more extensive numerical investigation. For the sake of com- 
pleteness we also show in figure [TO] similar data for the case 
where fcg'^ = 0. We find that qualitatively the data has a very 
similar behaviour to the case where kf^ — /c™. 




FIG. 9: (colour online) The two-time auto-correlation function C (r) 
for fci = 1, ^2 = 100. The insets shows the same data in a semi-log 
scale. Top: ^4 is varied with constant kg = 0.001. BOTTOM: fca is 
varied with constant = 0.1. 




FIG. 10: (colour online) The correlation function C(r) when fej^" = 
0. As previously, ki — 1 and k2 = 100. The data labelled as a 
corresponds to fc™ = 0.001 while b corresponds to k°^ = 0.01. For 
each a and b the superscripts 1, 2 and 3 refer to different values of 
k4 = 0.1, 1 and 10 respectively. The inset shows the same plot on a 
semi-log scale. 



B. Multiple coupled switches 

Many bacterial genomes contain multiple phase-varying 
genetic switches, which may demonstrate correlated flip- 
ping. For example, in uropathogenic E. coli, the Jim and pap 
switches, which control the production of different types of 
fimbriae, have been shown to be coupled Il42ll43l . Although 
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these two switches operate by different mechanisms, it is also 
likely that multiple copies of the same switch are often present 
in a single cell. This may be a consequence of DNA replica- 
tion before cell division (in fast-growing E. coli cells, division 
may proceed faster than DNA replication, resulting in up to 
^ 8 copies per cell). Randomly occurring gene duplication 
events, which are believed to be an important evolutionary 
mechanism, might also result in multiple copies of a given 
switch on the chromosome. It is therefore important to under- 
stand how multiple copies of the same switch would be likely 
to affect each other's function |44|. 

Let us suppose that there are two copies of our model switch 
in the same cell. Each copy contributes to and is influenced by 
a common pool of molecules of enzyme R. Our model is still 
described by the set of reactions ([T]l, but now the copy number 
of Son and S'off can vary between and 2 (with the constraint 
that the total number of switches is 2). 

To measure correlations between the states of the two 
switches (denoted si and S2) we define the two switch joint 
probability p2{si,t; S2,t') as the probability that switch 1 is 
in state si at time t and switch 2 is in state S2 at time t' . This 
function is the natural extension of the previously defined two- 



time probability for a single switch. Thus, in analogy to ( 60 1, 
we can define a two-time correlation function: 

P2(on, <;on, i + r) p2(off, off, i + r) 



C2(r) = 

Pon PoS 

(61) 

where pon (Pos) is again the steady-state probability for a 
single switch to be on (off). If the two switches are com- 
pletely uncorrected, we expect that p2{on,t;on,t') — p^^ 
and p2(off, off, t') — p^jj, so that C2(t) = (given that 
Pon + Poff — !)• In contrast, if the switches are completely 
con-elated, p2(on, on, i') = pon, _P2(off, <; off, <') = poff 
and C2(t) ~ 1. For completely anti-correlated switches, 
we expect that p2 (on, i; on, t') = p2(off, t; off, t') = 0, and 
C2(r) ~ 



-1. In Fig. 
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we plot the function C2(t) for two 
identical coupled switches, for several parameter sets. Our re- 
sults show that for small values of ^4, there is correlation be- 
tween the two switches, over a time period w lOk^^, which 
is of the same order as the typical time spent in the on state 
for these parameter values. Our results also show that the na- 
ture of these correlations depends strongly on kf^. In the case 
where kf^ 



11 



one can see that 

the correlation is positive, meaning that the two switches are 



fc™ (top panel of Fig. 



more likely to be in the same state 
set to zero (bottom panel of Fig 



In contrast, when kf^ is 
111, the correlation is neg- 



ative, meaning that the two switches are more likely to be in 
different states. 

To understand these correlations, consider the extreme sit- 
uation where both the two switches are off, and the number 
molecules of R has dropped to zero. In this case, the only pos- 
sible event is a fc4 mediated switching which could take place, 
for instance, for the first switch. Then, once the first switch 
is on, it will start producing more enzyme, and, if kf^ ^ 0, 
this will enhance the probability for the second switch to flip 
on too. This might explain why, when kf^ = fc™' we see a 
positive correlation between the two switches. On the other 
hand, if we consider the opposite situation where both the two 
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FIG. 1 1 : (colour online) Normalised two-time correlation function 
C2(t) for two identical switches. The parameter values are: ki — 1, 
k2 = 100,fc;^" = 0.001. In the top panel kf = fc^" while in the 
bottom panel kf^ = 0. The parameter k4 is varied from 0.1 to 100 
in each case. 



switches are on, and the number of molecules of R is around 
its plateau value, then the on to off switching probability for 
the two switches will be at its maximum. However, after one 
of the switches has flipped (e.g. the first), the switching prob- 
ability will start decreasing, this reducing the flipping rate for 
the second switch. This suggests that fc™ may have the effect 
of inducing negative correlations, while k'^^ induces positive 
correlations. We also point out the presence of a small peak 
in C2(t) in Fig. 



1 1 (indicated by the arrow) which suggests 



the presence of a time delay: when one switch flips, the other 
tends to follow a short time later We leave the detailed prop- 
erties of these correlations and their parameter dependence to 
future work. 



VI. SUMMARY AND OUTLOOK 

In this paper we have made a detailed study of a generic 
model of a binary genetic switch with linear feedback. The 
model system was defined in section II by the system of chem- 
ical reactions ([T]). Linear feedback arises in this switch be- 
cause the flipping enzyme R is produced only when the switch 
is in the on state, and the rate of flipping to the off state in- 
creases linearly with the amount of R. Thus, when the switch 
is in the on state the system dynamics inexorably leads to a flip 
to the off state. We have shown that this effect can produce a 
peaked flip time distribution and a bimodal probability distri- 
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bution for the copy number of R. A mean field description 
does not reproduce this phenomenology and so a stochastic 
analysis is required. 

We have studied this model analytically, obtaining exact 
solutions for the steady state distribution of the number of 
R molecules, as well as for the flip time distributions in the 
two different measurement ensembles defined in Section |IV] 
the Switch Change Ensemble and the Steady State Ensemble. 
We have shown how these ensembles are related and demon- 
strated that the flip time distribution in the Switch Change En- 
semble may exhibit a peak but the flip time distribution in the 
Steady State Ensemble can never do so. We also provide a 
generic relationship between the flip time distribution sam- 
pled in the two different ensembles. Given that in single-cell 
experiments, measuring the flip time distribution in the SCE 
is much more demanding than in the SSE, our result provides 
a way to access the SCE flip time distribution by making mea- 
surements only in the SSE. Our flip time calculations are rem- 
iniscent of persistence problems in non-equilibrium statistical 
physics where, for example, one is interested in the time an 
Ising spin stays in one state before flipping. However, be- 
cause of the linear feedback of our model switch, the flip time 
distribution is not expected to have a long tail as in usual per- 
sistence problems, rather it is the shape of the peak of the 
distribution which is of interest. 

By studying numerically the time correlations of a single 



switch, using the two time autocorrelator ( 60 1, we have shown 



that our model switch can play the role of a primitive "mem- 
ory module". The two time autocorrelator displays nontrivial 
behaviour including rather slow decay, which would be wor- 
thy of further study. We have also investigated the behaviour 
of two coupled switches within the same cell, and showed 
that both positive and negative correlations could be produced 
by choosing the parameters appropriately. In particular for 
fcg** = 0, as is the case for the Jim switch, anti-correlations 
were observed, implying that if one switch were on at time t, 
the other would tend to be off at that time and for a subsequent 
time of about one switch period. 

Many open questions and problems remain. At a technical 
level one would like to compute correlations of a single 
switch analytically and be able to treat the multiple switch 
system. The model itself could be refined in several ways, for 
example, by introducing nonlinear feedback fl31 l46l . It has 
been shown that such feedback allows nontrivial behaviour 
even at the level of a piecewise deterministic Markov process 
approximation [46!], where one assumes a deterministic 
evolution for the enzyme concentration, but a stochastic 
description for the switching. At present our model includes 
no explicit coupling to the environment, but such coupling 
could be included in a simple way by adding into the model 
environmental control of parameters or fc4. To make a 
closer connection to real biological switches, such as Jim, one 
could extend the model to include, for example, multiple and 
cooperative binding of the enzymes ll26l l27l . One particularly 
exciting direction, which we plan to pursue in future work, 
is to develop models for growing populations of switching 
cells, in which cell growth is coupled to the switch state. 
Such models could lead to a better understanding of the role 



of phase variation in allowing cells to survive and proliferate 
in fluctuating environments. 
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APPENDIX A: SOLUTION FOR THE SURVIVAL 
PROBABILITY 



We show here how to solve Eq.(44ai using the method of 



characteristics (see e.g. 1411 ). Introducing the new variable 

r{z, t), we set 

dhonizir),t{r)) _ dt 8 ~ ,..,dzd_~ 

dr ~ 9r9t ™^ ' ^+ 9r9z ' ^ 

= ^L,{z,t) + [h{z-l) + kl"z]—h,,{z,t) . (Al) 
ot az 

We can then identify the derivatives of t and z with respect to 
r as: 

dt dz , ^ 

- = 1, - = fci(z-l) + fc™z. (A2) 

Next, we solve these equations for t{r) and z{r) using initial 
conditions t(0) ~ and 2:(0) — zq: 

t{r) = r , z{r) = fciTo„ + e''^'^""{zo ~ fciTo„) , (A3) 

where Ton — {ki + fc™)^^. The reduced ordinary differential 
equation (ODE) for ho„ is: 



dr 



^ [k2{z{r) ~ 1) ~ kT]K,{r) , (A4) 



Substituting in the above relation z{r) with its expression 
given in ( |A3| l, we get an ordinary differential equation for 
hon{f), which can be solved by separation of variables: 



dhn 



{-k2kf - kT/T,, + e'-/^"fc2(2o/ro„ - k,)) dr . 



(A5) 

Solving the above equation using the initial condition hon{r = 
0) = W{zq), we arrive at 



honir) = exp 



W{zo) , (A6) 



where Won = fc™ + ^2(1 — ^iTon)- Substituting then from 

(A3 \ r ^ t and zq — > kiT^a + e~*/'^"(z — fciTon) one finally 



recovers d46 



16 



APPENDIX B: BACKWARDS EVOLUTION EQUATIONS 
FOR FLIP TIME DISTRIBUTION 



In this appendix we show how the result ( 55 i can be ob- 
tained by considering the backward survival probability: 



hs ino.t) = hs{n,Q\nQ, ~t) , 



(Bl) 



which is the probability that the system has survived in the 
state s without flipping and with n enzymes at time know- 
ing that it had riQ enzyme molecules at a past time —t. The 
probability h~ will verify the backward master equation 



d_ 
dt 



{no, t) = nnkih^ {no - l,t) + k^h^ {uo + 1, 

- (nofci + + uokl + kl)h-{no,t) . (B2) 



In section |IV] we used the forward master equation to com- 
pute the flip time distribution in two steps. First, we computed 
the forward survival probability hs{n,t) with two possible ini- 
tial conditions, to distinguish the two possible scenarios of 
measurement. Second, we summed this survival probability 
over all possible final configurations, and took the time deriva- 
tive in order to enforce a flipping at the end of the sampling. 

An analogous calculation (which we do not detail) can be 



carried out considering the backward master equation ( B2 1, 
and the final result has to be the same. In fact, we can con- 



sider the rh.s. of (B2 1 as a generator of the backward dynam- 
ics. Thus the solution of the backward evolution equation will 
have as boundary condition the statistics of the final config- 
uration at time 0, and will yield the statistics of the possible 
corresponding initial configurations at —t (with the additional 
constraint that the switch never flipped). Since for both SCE 
and SSE we condition that on switch flips at t = 0, the bound- 



ary condition of (B2 1 has to be taken when the switch is flip- 
ping from state s to state 1 — s, and thus corresponds to: 



h-{n,0) = W!^'',{n) 



(B3) 



step described above. The advantage is that now our boundary 
condition is the same for both the SCE and the SSE. 

We can relate hj to Fg by noting that hj {no, t) is the 
probability that the switch has not flipped going backward for 
a time t. We now have to made a distinction between the SCE 
and the SSE, since what happens at time —t is precisely the 
initial ensemble. For the case of the SCE, we want the switch 
to flip at time —t, therefore the flip time distribution is given 
by: 



SCE 



{T)^^dTY.h- 



{no,T) 



(B4) 



On the other hand, for the case of the SSE, there is no flipping 
at — t to enforce and the flip time distribution F|se simply 
proportional to the survival probability: 



(B5) 



irdT'Enohs{no,T') ' 
The denominator in ( |B5| l is chosen to ensure normalisation 

/dTFfE(r) = 1. 

Furthermore, we can compute the average flip time in the 



SCE using (B4 1 



dT' T'F^^^{T') 

= / dT' V/i;(no,T') , (B6) 

where an integration by parts has been performed. We can see 



then that the denominator in Eq.(B5 



flip time. Finally, integrating Eq.(B4 



is exactly the average 
from T to infinity and 



replacing the result in ( |B6| l, we obtain 

j^F!''^{r) dT' 



where W^^"--^ is defined in (40 1. This is the analogue of the first and the result 



/o°°T'F|CE(T') dT' ■ 
is recovered. 



(B7) 
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